
# ==============================================================================
# DESCRIPTION ----
# ==============================================================================

# Conducts robustness checks related to pre-trends for time-series cross
# sectional matching analysis of amakudari hires on government loans received by 
# private sector firms. 

# ==============================================================================
# LIBRARIES AND IMPORT ----
# ==============================================================================

# Options
options(scipen=999)

# Libraries
library(tidyverse)
library(DeclareDesign)
library(PanelMatch)

# Functions
source("code/0.2_functions.R")

# Data import
load("data/loans_amakudari.Rdata")

# ==============================================================================
# FINAL DATA PREP ----
# ==============================================================================

# Some duplicates in data
financials <- loans_covs %>% 
  distinct(firm_dest_en, year, .keep_all = TRUE) %>%
  filter(year != 2020)

# Create alternative amakudari coding where 1 for all periods after hire
fin <- financials %>% 
  mutate(amakudari_fill = na_if(as.numeric(amakudari_binary), 0),
         amakudari_match = na_if(as.numeric(amakudari_binary), 0)) %>%
  group_by(nikkei_code) %>%
  fill(amakudari_fill, .direction = "down") %>%
  fill(amakudari_match, .direction = "downup") %>%
  ungroup() %>%
  mutate(amakudari_fill = as.numeric(replace_na(amakudari_fill, 0)),
         amakudari_match = as.numeric(replace_na(amakudari_match, 0)),
         year = as.integer(as.numeric((as.factor(year)))),
         nikkei_code = as.numeric(factor(nikkei_code)))

fin <- as.data.frame(fin)

# Create amakudari codings for each ministry -----------------------------------
ministry <- financials %>% 
  mutate(
    amakudari_binary_METI_MOF = amakudari_binary_METI + amakudari_binary_MOF,
    amakudari_binary_METI_MOF = ifelse(amakudari_binary_METI_MOF >= 1, 1, amakudari_binary_METI_MOF),
    amakudari_binary_other = amakudari_binary_MIAC + amakudari_binary_MLIT + 
      amakudari_binary_MHLW + amakudari_binary_MAFF + amakudari_binary_MEXT +
      amakudari_binary_MOD + amakudari_binary_MOE + amakudari_binary_MOFA +
      amakudari_binary_MOJ,
    amakudari_binary_other = ifelse(amakudari_binary_other >= 1, 1, amakudari_binary_other),
    across(starts_with("amakudari"), ~na_if(as.numeric(.), 0))
  ) %>%
  group_by(nikkei_code) %>%
  fill(c(amakudari_n_Other:amakudari_binary_other), .direction = "down") %>%
  ungroup() %>%
  mutate(
    across(starts_with("amakudari"), ~as.numeric(replace_na(., 0))),
    year = as.integer(as.numeric((as.factor(year)))),
    nikkei_code = as.numeric(factor(nikkei_code))
  )
ministry <- as.data.frame(ministry)

# ==============================================================================
# RESTRICTING MATCHES TO ADDITIONAL YEARS AND PLACEBO TESTS: ALL MINISTRIES ----
# ==============================================================================
# 2 lag estimates --------------------------------------------------------------
loan_result <- run_panelmatch(
  data = fin, treatment = "amakudari_fill", outcome = "debt_total_public", 
  lag = 2)

# FIGURE A14
ggplot(loan_result$cleaned_result_table) + gglayers_update +
  scale_y_continuous(limits = c(-5000, 20000), breaks= seq(-5000, 20000, 5000))
ggsave("figures/a14_tscs_loan_2lags.pdf", width = 10, height = 5)

# TABLE A11: Output table
stargazer::stargazer(
  loan_result$cleaned_result_table, 
  summary = FALSE, digits = 2, rownames = FALSE,
  title = "Estimated effect of bureaucratic hires on size of government loans received, requiring 2 lag periods",
  label = "tab: tscs_loan_lag",
  out = "tables/a11_tscs_loan_lag.tex",
  notes = paste("Note: Matched sets =", length(loan_result$estimate_object$matched.sets))
)

# Placebo test -----------------------------------------------------------------
# Open a PDF device
pdf("figures/a15_tscs_loan_2lags_placebo.pdf", width = 10, height = 5) 

# Run placebo test for lag years
placebo_test(
  pm.obj = loan_result$match_object, panel.data = loan_result$panel_object, 
  number.iterations = 1000,
  plot = TRUE, se.method = "bootstrap", parallel = T, num.cores = 4)

# Close the PDF device
dev.off()

# ==============================================================================
# RESTRICTING MATCHES TO ADDITIONAL YEARS AND PLACEBO TESTS: METI & MOF ----
# ==============================================================================
# Requiring lags ---------------------------------------------------------------
# 2 lags
loan_result_meti_mof2 <- run_panelmatch(
  data = ministry, treatment = "amakudari_binary_METI_MOF", 
  outcome = "debt_total_public", lag = 2)

# 3 lags
loan_result_meti_mof3 <- run_panelmatch(
  data = ministry, treatment = "amakudari_binary_METI_MOF", 
  outcome = "debt_total_public", lag = 3)

# 4 lags
loan_result_meti_mof4 <- run_panelmatch(
  data = ministry, treatment = "amakudari_binary_METI_MOF", 
  outcome = "debt_total_public", lag = 4)

# Plot results -----------------------------------------------------------------
# FIGURE A16b: 2 lags
ggplot(loan_result_meti_mof2$cleaned_result_table) + gglayers_update
ggsave("figures/a16b_tscs_loan_meti_mof_2.pdf", width = 10, height = 5)

# FIGURE A16c: 3 lags
ggplot(loan_result_meti_mof3$cleaned_result_table) + gglayers_update
ggsave("figures/a16c_tscs_loan_meti_mof_3.pdf", width = 10, height = 5)

# FIGURE A16d: 4 lags
ggplot(loan_result_meti_mof4$cleaned_result_table) + gglayers_update
ggsave("figures/a16d_tscs_loan_meti_mof_4.pdf", width = 10, height = 5)

# Create table for 2 lag periods -----------------------------------------------
# TABLE A12
stargazer::stargazer(
  loan_result_meti_mof2$cleaned_result_table, 
  summary = FALSE, digits = 2, rownames = FALSE,
  title = "Estimated effect of bureaucratic hires on size of government loans received, METI and MOF re-hires only, requiring 2 lag periods",
  label = "tab: tscs_loan_lag",
  out = "tables/a12_tscs_loan_lag_meti_mof.tex",
  notes = paste("Note: Matched sets =", length(loan_result$estimate_object$matched.sets))
)

# Run placebo tests for lag years ----------------------------------------------
# FIGURE A17a: 2 lag years
pdf("figures/a17a_tscs_loan_meti_mof_2_placebo.pdf", width = 10, height = 5) 
placebo_test(
  pm.obj = loan_result_meti_mof2$match_object, 
  panel.data = loan_result_meti_mof2$panel_object, number.iterations = 1000,
  plot = TRUE, se.method = "bootstrap", parallel = T, num.cores = 4)
dev.off()

# FIGURE A17b: 3 lag years
pdf("figures/a17b_tscs_loan_meti_mof_3_placebo.pdf", width = 10, height = 5) 
placebo_test(
  pm.obj = loan_result_meti_mof3$match_object, 
  panel.data = loan_result_meti_mof3$panel_object, number.iterations = 1000,
  plot = TRUE, se.method = "bootstrap", parallel = T, num.cores = 4)
dev.off()

# FIGURE A17c: 4 lag years
pdf("figures/a17c_tscs_loan_meti_mof_4_placebo.pdf", width = 10, height = 5) 
placebo_test(
  pm.obj = loan_result_meti_mof4$match_object, 
  panel.data = loan_result_meti_mof4$panel_object, number.iterations = 1000,
  plot = TRUE, se.method = "bootstrap", parallel = T, num.cores = 4)
dev.off()

# ==============================================================================
# PRE-TREATMENT TRENDS: ALL BUREAUCRATS ----
# ==============================================================================
# One lag period matching for all bureaucratic appointments 
loan_result_pre <- run_panelmatch(
  data = fin, treatment = "amakudari_fill", outcome = "debt_total_public")

# Treated firms ----------------------------------------------------------------
# Pull treated firms out of matched set
loan_match.msets <- extract(loan_result_pre$match_object)
treated_firms_years <- names(loan_match.msets)
treated_firms <- sub("\\..*", "", treated_firms_years)

# Get the first year of amakudari hire for each firm
year_hire <- as.integer(sub(".*\\.(.*)", "\\1", treated_firms_years))
treated_firms_years <- as.data.frame(cbind(treated_firms, year_hire)) %>%
  rename(nikkei_code = treated_firms)

# Join back and calculate relative year and mean
normalized_data <- fin %>%
  filter(nikkei_code %in% treated_firms, amakudari_fill == 0) %>%
  mutate(nikkei_code = as.character(nikkei_code)) %>%
  left_join(treated_firms_years, by = "nikkei_code") %>%
  mutate(relative_year = year - as.integer(year_hire))

# Aggregate across firms by relative year
relative_year_mean <- normalized_data %>%
  group_by(relative_year) %>%
  summarize(
    `Government loans` = mean(debt_total_public), 
    `Assets` = mean(total_assets, na.rm = T), 
    `Liabilities` = mean(total_liabilities, na.rm = T), 
    `Employees` = mean(employees, na.rm = T), 
    `EBITDA` = mean(ebitda, na.rm = T), 
    `Revenue` = mean(operating_revenue, na.rm = T), 
    `Profit` = mean(gross_profit, na.rm = T), 
    .groups = "drop"
  )

# Weighted mean for controls ---------------------------------------------------
# Extract matched set metadata
loan_match.msets <- extract(loan_result_pre$match_object)
treated_firms_years <- names(loan_match.msets)
treated_firms <- sub("\\..*", "", treated_firms_years)
year_hire <- as.integer(sub(".*\\.(.*)", "\\1", treated_firms_years))

# Build a data frame of matched control units
matched_controls_df <- purrr::map2_dfr(
  .x = loan_match.msets,
  .y = treated_firms_years,
  .f = function(control_ids, treated_key) {
    treated_firm <- sub("\\..*", "", treated_key)
    year_hire <- as.integer(sub(".*\\.", "", treated_key))
    tibble(
      nikkei_code = as.character(control_ids),
      treated_firm = treated_firm,
      year_hire = year_hire,
      treated_key = treated_key
    )
  }
)

# Extract weights
weights_df <- weights(loan_match.msets) %>%
  enframe(name = "treated_key", value = "weights_list") %>%
  unnest_longer(weights_list, indices_to = "nikkei_code", values_to = "weight") %>%
  mutate(nikkei_code = as.character(nikkei_code))

# Merge weights into matched control df
matched_controls_df <- matched_controls_df %>%
  left_join(weights_df, by = c("treated_key", "nikkei_code"))

# Join to panel data and compute relative year
control_normalized_data <- fin %>%
  filter(amakudari_fill == 0) %>%
  mutate(nikkei_code = as.character(nikkei_code)) %>%
  inner_join(matched_controls_df, by = "nikkei_code") %>%
  mutate(relative_year = year - as.integer(year_hire))

# Weighted aggregation by relative year
control_relative_year_weighted_mean <- control_normalized_data %>%
  group_by(relative_year) %>%
  summarize(
    `Government loans` = weighted.mean(debt_total_public, weight, na.rm = TRUE),
    `Assets` = weighted.mean(total_assets, weight, na.rm = TRUE),
    `Liabilities` = weighted.mean(total_liabilities, weight, na.rm = TRUE),
    `Employees` = weighted.mean(employees, weight, na.rm = TRUE),
    `EBITDA` = weighted.mean(ebitda, weight, na.rm = TRUE),
    `Revenue` = weighted.mean(operating_revenue, weight, na.rm = TRUE),
    `Profit` = weighted.mean(gross_profit, weight, na.rm = TRUE),
    .groups = "drop"
  )

# Create plots of means --------------------------------------------------------
# Add group labels to each data frame
treated_df <- relative_year_mean %>%
  mutate(group = "Treated")

control_weighted_df <- control_relative_year_weighted_mean %>%
  mutate(group = "Control")

# Reshape to long format for plotting
combined_df <- bind_rows(treated_df, control_weighted_df) %>%
  pivot_longer(cols = -c(relative_year, group), 
               names_to = "variable", values_to = "value") %>%
  filter(relative_year < 0) %>%
  mutate(
    group = factor(group, levels = c("Treated", "Control")),
    variable = factor(variable, 
                      levels = c("Government loans", "Assets", "Liabilities",
                                 "Employees", "EBITDA", "Revenue", "Profit"))
    )

# FIGURE A25
ggplot(combined_df, aes(x = relative_year, y = value, linetype = group)) +
  geom_line(size = 0.5) +
  facet_wrap(~ variable, scales = "free_y", ncol = 2) +
  labs(x = "Relative year before treatment", 
       y = "Mean Value", color = "black", linetype = NULL) +
  scale_linetype_manual(values = c("Treated" = "solid", "Control" = "dashed")) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "bottom",
    strip.text = element_text(size = 13, face = "bold")
  )

ggsave("figures/a25_pre_trends.pdf", width = 10, height = 7)

# ==============================================================================
# PRE-TREATMENT TRENDS: METI/MOF ----
# ==============================================================================
# One lag period matching for all meti/mof appointments 
loan_result_pre_meti_mof <- run_panelmatch(data = ministry, 
    treatment = "amakudari_binary_METI_MOF", outcome = "debt_total_public")

# Treated firms ----------------------------------------------------------------
# Pull treated firms out of matched set
loan_match.msets <- extract(loan_result_pre_meti_mof$match_object)
treated_firms_years <- names(loan_match.msets)
treated_firms <- sub("\\..*", "", treated_firms_years)

# Get the first year of amakudari hire for each firm
year_hire <- as.integer(sub(".*\\.(.*)", "\\1", treated_firms_years))
treated_firms_years <- as.data.frame(cbind(treated_firms, year_hire)) %>%
  rename(nikkei_code = treated_firms)

# Join back and calculate relative year and mean
normalized_data <- fin %>%
  filter(nikkei_code %in% treated_firms, amakudari_fill == 0) %>%
  mutate(nikkei_code = as.character(nikkei_code)) %>%
  left_join(treated_firms_years, by = "nikkei_code") %>%
  mutate(relative_year = year - as.integer(year_hire))

# Aggregate across firms by relative year
relative_year_mean <- normalized_data %>%
  group_by(relative_year) %>%
  summarize(
    `Government loans` = mean(debt_total_public), 
    `Assets` = mean(total_assets, na.rm = T), 
    `Liabilities` = mean(total_liabilities, na.rm = T), 
    `Employees` = mean(employees, na.rm = T), 
    `EBITDA` = mean(ebitda, na.rm = T), 
    `Revenue` = mean(operating_revenue, na.rm = T), 
    `Profit` = mean(gross_profit, na.rm = T), 
    .groups = "drop"
  )

# Weighted mean for controls ---------------------------------------------------
# Extract matched set metadata
loan_match.msets <- extract(loan_result_pre_meti_mof$match_object)
treated_firms_years <- names(loan_match.msets)
treated_firms <- sub("\\..*", "", treated_firms_years)
year_hire <- as.integer(sub(".*\\.(.*)", "\\1", treated_firms_years))

# Build a data frame of matched control units
matched_controls_df <- purrr::map2_dfr(
  .x = loan_match.msets,
  .y = treated_firms_years,
  .f = function(control_ids, treated_key) {
    treated_firm <- sub("\\..*", "", treated_key)
    year_hire <- as.integer(sub(".*\\.", "", treated_key))
    tibble(
      nikkei_code = as.character(control_ids),
      treated_firm = treated_firm,
      year_hire = year_hire,
      treated_key = treated_key
    )
  }
)

# Extract weights
weights_df <- weights(loan_match.msets) %>%
  enframe(name = "treated_key", value = "weights_list") %>%
  unnest_longer(weights_list, indices_to = "nikkei_code", values_to = "weight") %>%
  mutate(nikkei_code = as.character(nikkei_code))

# Merge weights into matched control df
matched_controls_df <- matched_controls_df %>%
  left_join(weights_df, by = c("treated_key", "nikkei_code"))

# Join to panel data and compute relative year
control_normalized_data <- fin %>%
  filter(amakudari_fill == 0) %>%
  mutate(nikkei_code = as.character(nikkei_code)) %>%
  inner_join(matched_controls_df, by = "nikkei_code") %>%
  mutate(relative_year = year - as.integer(year_hire))

# Weighted aggregation by relative year
control_relative_year_weighted_mean <- control_normalized_data %>%
  group_by(relative_year) %>%
  summarize(
    `Government loans` = weighted.mean(debt_total_public, weight, na.rm = TRUE),
    `Assets` = weighted.mean(total_assets, weight, na.rm = TRUE),
    `Liabilities` = weighted.mean(total_liabilities, weight, na.rm = TRUE),
    `Employees` = weighted.mean(employees, weight, na.rm = TRUE),
    `EBITDA` = weighted.mean(ebitda, weight, na.rm = TRUE),
    `Revenue` = weighted.mean(operating_revenue, weight, na.rm = TRUE),
    `Profit` = weighted.mean(gross_profit, weight, na.rm = TRUE),
    .groups = "drop"
  )

# Create plots of means --------------------------------------------------------
# Add group labels to each data frame
treated_df <- relative_year_mean %>%
  mutate(group = "Treated")

control_weighted_df <- control_relative_year_weighted_mean %>%
  mutate(group = "Control")

# Reshape to long format for plotting
combined_df <- bind_rows(treated_df, control_weighted_df) %>%
  pivot_longer(cols = -c(relative_year, group), 
               names_to = "variable", values_to = "value") %>%
  filter(relative_year < 0) %>%
  mutate(
    group = factor(group, levels = c("Treated", "Control")),
    variable = factor(variable, 
                      levels = c("Government loans", "Assets", "Liabilities",
                                 "Employees", "EBITDA", "Revenue", "Profit"))
  )

# FIGURE A26
ggplot(combined_df, aes(x = relative_year, y = value, linetype = group)) +
  geom_line(size = 0.5) +
  facet_wrap(~ variable, scales = "free_y", ncol = 2) +
  labs(x = "Relative year before treatment", 
       y = "Mean Value", color = "black", linetype = NULL) +
  scale_linetype_manual(values = c("Treated" = "solid", "Control" = "dashed")) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "bottom",
    strip.text = element_text(size = 13, face = "bold")
  )

ggsave("figures/a26_pre_trends_meti_mof.pdf", width = 10, height = 7)
